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Abstract - The dynamics of a confined fluid of Bose atoms is treated within the linear response 
regime, with a view to establishing a current-density functional formalism for an 

On 

i— h inhomogeneous superfluid state. After evaluating in full detail a simplified case of an external 

coupling to the density and phase of the condensate, the theory is extended to include the 

coupling to the total current density. The Kohn-Sham response functions of the condensate and 

all the exchange-correlation kernels for the superfluid are introduced from the microscopic 

^ equations of motion and are expressed in a physically transparent way through functional 
O 

00 derivatives of correlation functions. A microscopic formula for the superfluid density is derived 
G\ 

+^ and used to introduce a generalized hydrodynamic approach for a weakly inhomogeneous two- 

fluid model in isothermal conditions. Local-density expressions are thereby derived for the 

£h velocities of first and second sound in the weakly inhomogeneous superfluid and for visco- 

O 

elastic functions describing the transition from the hydrodynamic to the collisionless regime. 
Landau's hydrodynamic theory and known results in Green's functions language are recovered 
in the limiting case of a homogeneous superfluid. 
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1. Introduction 

A considerable experimental effort is currently being devoted to the elucidation of the 
dynamical behaviour of dilute alkali vapours which have undergone Bose-Einstein condensation 
in magnetic traps. These experimental studies have concerned the excitation of low-lying shape- 
deformation modes [1 - 3], the propagation of sound waves in the condensate and its thermal 
cloud [4, 5] and antisymmetric oscillations of the condensate and the thermal cloud 
corresponding to second sound [5]. The experimental progress has stimulated a number of 
theoretical approaches involving quantal hydrodynamic descriptions of confined Bose- 
condensed fluids [6 - 11]. The relationship between the dynamics of such inhomogeneous 
fluids and the well-known dynamic behaviour of a homogeneous Bose superfluid [12 - 14] also 
is of considerable interest. 

The developments that we have recalled above have motivated the present theoretical 
study. We bring forward the proposal made in our earlier work [7] to extend to a superfluid of 
neutral Bose particles the recent progress made in developing a time-dependent current 
functional theory for the dynamics of inhomogeneous electron systems. For a many-electron 
system in the normal state the basis of the theory comes from a set of theorems by Runge and 
Gross [15], showing that the problem of its dynamical behaviour in an external time-dependent 
potential can be mapped into that of non-interacting electrons in a self-consistent effective 
potential. These general theorems have been extended to inhomogeneous superconductors [16], 
the central result being a set of time-dependent Bogolubov - de Gennes equations which include 
exchange and correlation (xc) effects. In actual applications to electronic systems in the normal 
state (for reviews see Gross and Kohn[17] and Gross et al. [18]) low-frequency phenomena 
have been usefully described by means of the so-called "adiabatic local density approximation" 
(ALDA), in which the dynamic xc potential is evaluated as in the corresponding static problem 
from the xc energy density of the homogeneous electron gas at the local time-dependent density. 
More recently the search for a fully dynamic local-density approximation has led Vignale and 
Kohn [19, 20] to develop a current-density formulation of the theory. In brief, these authors 
have shown that the linear response of a weakly inhomogeneous electron system to an external 
time-dependent potential can be expressed in terms of (i) a dynamic xc vector potential built 
from xc kernels which are taken from the homogeneous fluid at the local equilibrium density, 
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and (ii) a Kohn-Sham current response matrix to be evaluated from a suitable set of single- 
particle orbitals. Such a framework supplements the ALDA in embodying not only the electron- 
gas compressibility but also plasmon dispersion and damping as well as transverse-current 
fluctuations [21, 22] and allows a unified treatment of the damping of collective excitations 
from the Landau and mode-coupling mechanisms. The formalism can be readily reduced to 
yield generalized hydrodynamic equations showing that the xc kernels have the meaning of 
frequency-dependent viscosity coefficients and elastic moduli [23]. 

In the case of a normal fluid at zero temperature as treated by Vignale and Kohn [19] the 
only relevant dynamical variable is the current density. For an extension of this approach to a 
superfluid of neutral Bose particles in isothermal conditions at finite temperature it is 
immediately realized that the set of basic dynamical variables must include, in addition to the 
total current density, the density of the condensate and its phase, the gradient of the latter giving 
the superfluid velocity field below threshold for vortex generation. In section 2 we present in 
full detail the theory of the linear response of the two dynamical variables of an inhomogeneous 
condensate to a dynamic gauge-breaking external field, starting from the microscopic equation 
of motion for the order parameter in terms of the condensate self-energy. This response is 
explicitly shown to have the Hohenberg-Kohn-Sham structure compatible with a mapping of 
the interacting system into a single-particle reference system. We thus identify the Kohn-Sham 
response functions of the condensate as well as the relevant xc kernels, the latter being 
expressed in the form of functional derivatives of correlation functions. The connection between 
the condensate kernels and the Green's functions which are more commonly used to describe 
the dynamics of homogeneous superfluids [12, 24-27] is exposed in Appendix A. 

The response of the noncondensate and of the total current density, and the cross- 
couplings between condensate and noncondensate are then evaluated in section 3 and in 
Appendix B. A useful result which emerges from this treatment is a microscopic expression for 
the superfluid density. This introduces a two-fluid model for the generalized hydrodynamics of 
an inhomogeneous superfluid in the linear-response, weak inhomogeneity regime, which is 
developed in § 4. The crucial step is an extension of Landau's linearized hydrodynamic theory 
for the homogeneous superfluid [13] to finite-frequency phenomena, on the basis of the 
memory function formalism [28]. The basic assumption of such a theory is that both the 
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equilibrium density profile and the external perturbing fields are slowly varying in space, on 
length scales set by the interatomic distance and by c / ft) where c is the local speed of sound 
and ft) is the excitation frequency. Section 4 also presents a brief discussion of the evaluation of 
the xc kernels for a weakly inhomogeneous superfluid in a collisionless regime. Finally, section 
5 gives a summary of our main results. 

2. Response of the condensate 

We treat in this section the linear response of an inhomogeneous condensate to an 
external perturbation described by the Hamiltonian 



where y/(r) is the field operator and 7](r, t) is a symmetry -breaking scalar field. An account of 
the response of the noncondensate is deferred to section 3. 

2.1 Coupling to fluctuations in condensate density and superfluid velocity. 

As a first step in making explicit the physical meaning of the Hamiltonian (2.1), we 
follow Hohenberg and Martin [12] in introducing the transformation 



with h c (r) the density operator and <p(r) the phase operator of the condensate. Then, writing 
y/(r) =< y(r) > eq +Sy/(r) and with the definition n c (r) = \< y/(r) > eq \ 2 for the equilibrium 

density of the condensate, eqn (2.1) can be rewritten as 



where 5h c (r) = 2Re[< yr(r) > eq 8yr\r)] and S<p(r) = -Im[8yr\r)/ < i/A + (r) > eq ]. The fields 
a(r,t) = [n c (r)]" ] Re[< y/\r) > eq 7](r,t)] and &(r,t) = 2Im[< y/(r) > eq rj(r,t)] in eqn (2.3) 



evidently represent external couplings to fluctuations in the density and phase of the condensate. 
We use the symbol < ■ ■ ■ > to denote expectation values on the equilibrium ensemble at given 
temperature, with the suffix eq added to denote properties of the unperturbed fluid. 
The superfluid velocity operator v s (r) is introduced as 




(2.1) 



VA t (r) = [n e (r)] 1/2 exp[-/^(r)] 



(2.2) 




(2.3) 



v j (r) = m"V r (5^(r) 



(2.4) 



with m the particle mass [12]. We are assuming that the superfluid velocity vanishes 



everywhere at equilibrium. The perturbing Hamiltonian (2.3) is then rewritten as 



where 



H (t) = \d 3 r[a(r, t)5h c (r) + A(r, t) ■ v,(r)] , (2.5) 
A(r,0 = f L JjV r ^-V r ^(i J ,0 . (2.6) 
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In summary, the Hamiltonian (2.1) describes external fields acting on the density and 
phase fluctuations of an inhomogeneous condensate, the gradient of the phase fluctuations 
being related by eqn (2.4) to fluctuations in the superfluid velocity. 

2.2 Linear response functions and their connection to Green 's functions. 

A two-by-two matrix of (time-ordered) response functions for the condensate density 
and phase is introduced through the definitions 

z„w ) - d g^] > =< nm)m' )] > (2.7) 

etcetera. Here, 1 = (r^/J and the operators are in the Heisenberg representation. 

With the notation \jf(Y) = ¥l (l) and y/' (1) = y/' 2 (l) , we introduce symmetrized and 
equilibrium- weighed single-particle Green's functions G^(l,l' ) for the noncondensate as 

G^\U)=h- < ¥'M> eq G„/U')< y,(r r )>„ race] (2.8) 

where 

G ap (l,V ) = -i{< T[ ¥a (l) ¥ l(V )]>-< ¥a (l) >< vrj(l' ) >} . (2.9) 

The following expressions are then obtained for the time-ordered response functions from the 
definitions of 8h c (r) and <5<p(r): 

X„0,V ) = i[2n c ( ri )n c (r Y KWnW ) " G£>(U' )] , (2-10) 

XncHc (i,r ) = 2i[G 1 ^ ) (i,r ) + G{ 2 +) (i,i' )] , (2.1 1) 

X„Jl,V ) = -[n c (r v )TWu\W ) ~ G ( -\l,Y )] (2.12) 

and 

^ (1,1' ) = [n^ViG^&Y ) + G\~\U )] • (2.13) 
Of course, the expressions (2.10), (2.12) and (2.13) can be converted with the help of eqn 
(2.4) into expressions for response functions involving the superfluid velocity. 

In the present case the off-diagonal response functions coupling the superfluid velocity 
to the condensate density reflect the non-conservation of the number of particles in the 
condensate: by exciting fluctuations in the condensate one changes its density and hence induces 
a superfluid flow. Whereas in standard treatments of superfluid hydrodynamics the amplitude 
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fluctuations of the condensate can be neglected (see e. g. [14]), all the four response functions 
(2.10) - (2.13) are taken into account in our treatment, which is not restricted to low frequency. 

The symmetry properties of the response functions introduced above are easily assessed 
by a standard analysis of the spectral functions associated with the causal analogues of the 
Green's functions (2.8). Assuming a time-inversion-invariant unperturbed state, one can in 
particular prove the symmetry relation lmX fn (r,i? ;ft>) = -lmx n ,r,-(0), both functions 

being purely imaginary and even in 0) . 

2.3 Microscopic expressions for the response functions. 

We proceed to evaluate the response functions from the microscopic equation of motion 
for the expectation value < y/(l) > of the field operator in the presence of the perturbation 
H (t) . This is 

}J2G- 1 (1,2)< yr(2)>= 77(1) + cr(l) (2.14) 
where G is the free-particle Green's function, 

G- \1,V ) = V-^- + ^-V 2 Tl ~ V(r x ) + ) (2.15) 
ot l 2m 

with V (15 ) the confining potential and fi the chemical potential, and 

a(l) = J J2v( ri ,r 2 )< ^'(2)^(2)^(1) > (2.16) 

is the condensate self-energy, with v(r 1; r r ) the interparticle pair potential. 

By analogy with the definitions of the external fields entering eqn (2.3), we introduce 
the excess (ex) potentials 

a„(l) = [n c (r 1 )r 1 Re[< w'\^)> eq ofl)] (2-17) 

and 

^ x (l) = 2Im[<VA i (r 1 )> e?C 7(l)] . (2.18) 
The ex kernel a ex may be viewed as the shift in chemical potential due to the interactions in the 
fluid away from equilibrium. On the other hand, i} ex determines the violation of the continuity 
equation for the condensate through couplings to the noncondensate, as may be seen from the 
microscopic equation of motion for the condensate density fluctuations (see also Hohenberg and 
Martin [12]). 

Equation (2. 14) can then be written as 

< v(D>= JJ2G (1,2)< Y\%)> eq {[a(2) + a„(2)] 
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+ ^'[n c (r 2 )r 1 [^(2)+^,(2)]} . (2.19) 

This form of the equation of motion for the field operator can now be used to evaluate the 
response functions by taking functional derivatives with respect to the external fields in the 
linear response limit. 

Let us consider first the phase-phase response function ) . From its definition 

we have 



where 



+ G { ~\1,2) T } (2.20) 
8#(V ) 



G ( ± \l,V) = k< ys\r 1 )> eq G (l,V)<iff(r v )> eq ±c.c] . (2.21) 



It is evident from eqn (2.20) that G ( +) is related to an ideal phase-phase response function x^> 
aside from weighing factors determined by the equilibrium density of the real fluid: 

4; ) (i,r)=/[2n c (r 1 )n c (r 1 ,)r 1 G^ ) ar) . (2.22) 

A similar calculation for the phase-density response function attributes to Gq _) in eqn (2.21) the 

meaning of an ideal phase-density response function, 

4: ) (l,r) = [n c (r 1 )r 1 G<- ) (l,r) • (2.23) 

Equation (2.20) can therefore be written in its final form, 
y =yW+ry<°>® S ^«\ +y<°> ® 5<X ™ | l®y 

Aj <p(p Aj (p<p \-As(p<p ^* l n c ™<pn c ^ n c J A (pep 

+^®^K+XZ • (2-24) 

In eqn (2.24) we have omitted the arguments of all the functions involved and indicated with the 

symbol ® an integration over intermediate variables. Of course, eqn (2.4) immediately yields 
^ V! (l,r) = m- 2 V r V ri ^ w (l,r) . (2.25) 

The microscopic structure of the matrix of response functions is now clear. With the 
further definitions 

^(l,l') = -[n c (r 1 ,)r 1 G<- ) ai') (2.26) 

and 

^ e (l,r) = 2iG«(l,r) (2.27) 
we construct an ideal response matrix x (0) '• 



x (0) 



8 

fy(0) (0) \ 

A/n c n c A/n c <p 
y(0) y(0) 



(2.28) 



We also have the matrix K of excess kernels arising from the interactions: 
J5a„/8n e \ r 8a ex l8cp\\ 

Then the desired form of the matrix % of response functions in the real Bose fluid is given by 
X = X W +X W ®K®X (2-30) 

or 

X = U-X W ®KY l ®x W ■ (2-31) 
As a final remark we explicitly note that, at variance from what is customary for electron fluids, 
the kernels 8a ex I 8n c \ 9 and 8a ex I 8(p\ n in eqn (2.29) include Hartree contributions in addition 

to the exchange and correlation terms. 

The relevance of equation (2.30) to a density-functional treatment of the dynamics of an 
inhomogeneous superfluid will be discussed in section 2.4 immediately below. An alternative 
derivation of eqn (2.30), which highlights the meaning of the excess kernels (2.29) as single- 
particle self-energies for the real fluid and relates the present treatment to the work of Wong and 
Gould [25] on homogeneous Bose fluids, is given in Appendix A. 

2.4 Connection with time -dependent density functional theory 

The structure of the microscopic equations (2.30) for the response of the condensate at 
fixed noncondensate is a consequence of the fact that the self-energy o(l) is a function of only 

one space-time variable. In our derivation of these equations in § 2.3 we chose to use the 
response functions (2.28) as built from the ideal-gas Green's function G , in order to facilitate 

comparisons with earlier work on the homogeneous fluid. Hence, the form of eqn (2.30) is that 
appropriate for a Hohenberg-Kohn-Sham treatment, with a Kohn-Sham reference fluid which 
is, however, an ideal Bose gas except for weighing factors determined by the true equilibrium 
condensate density. 

A more appropriate set of Kohn-Sham response functions for a time-dependent density 
functional treatment of the condensate response can be obtained from a single-particle Green's 
function G KS (l,Y) defined through 

G-'(l,l') = G- 1 (l,r) *(U') • (2.32) 



< V(r,)> 



eq 
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The equilibrium self-energy term which has been included in eqn (2.32) needs to be subtracted 

from the RHS of eqn (2.14), which now reads 

f d2G K l(l,2) < y/(2) >= 77(1) + AoQ) (2.33) 

where 

4<7(1) = (7(1) - a eq {r x ) < > I < yr(r,) > eq (2.34) 

contains only dynamical effects. All other equations of § 2.3 remain formally the same, after 
replacement of G with G KS and of (7(1) with 4(7(1) . In particular, the matrix of Kohn-Sham 
response functions x KS 1S to be built from the Green's functions 

G%(l,V) = k< W'\^)> eq G KS {\,Y)<y{v v )> eq ±c.c] (2.35) 

and is in general non-diagonal. We also remark that a static mean-field potential due to the 
interactions has been included in earlier treatments of the dynamics of an inhomogeneous 
condensate in the random phase approximation [29, 30]. 

Thus, the Kohn-Sham reference system that we are proposing for a density functional 
approach to the dynamical xc effects in a superfluid explicitly contains the interactions through 
the equilibrium value of the condensate self-energy. Of course, expansion of the Kohn-Sham 
response functions in an appropriate basis set will be needed (see e.g. the suggestion of a 
Bogolubov - de Gennes basis set made by Wacker et al. [16] for a superconductor). On the 
other hand, we should emphasize that we are not advocating a density functional approach to 
the evaluation of the equilibrium state of the superfluid. In particular, a thermodynamic 
treatment based on the Gross-Pitaevskii equation is in quantitative agreement with experiment 
for confined vapours of alkali atoms which have undergone Bose-Einstein condensation [31]. 

We conclude this section by recalling that simple approximations on 4(7(1) , leading to a 
time-dependent Gross-Pitaevskii equation, have been found to be useful in the recent literature 
on collective excitations of the condensate (see e.g. [32]). 

3. Response of the noncondensate and of the whole fluid 

Having dealt with the linear response properties of the condensate in the preceding 

section, we now extend the treatment to the linear response of the whole Bose fluid. The 

perturbing Hamiltonian (2.1) is supplemented by the perturbation 

H"(t)= fd 3 rfd 3 r'ljf\r')U(r,r';t)li/(r) , (3.1) 
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involving a non-local scalar potential U(r,r' ;t) as needed to describe the current response. We 
shall first focus on the response of the noncondensate to this potential, before turning to treat 
the response of the whole fluid to the set of external perturbations given by ?7(r, t) and 
U(r,r';t). 

In earlier work on Bose fluids by the so-called dielectric formalism [12, 14, 24, 25], a 
diagrammatic analysis has often been used to divide the proper current response function into a 
'regular' part and a 'singular' part. The regular part is described by the set of irreducible 
diagrams and may be identified with the response of the normal component. The singular part 
can be expressed in terms of the response of the condensate, multiplied by a vertex function 
which contains off-diagonal (condensate-noncondensate) kernels. The same vertex function 
relates the condensate-noncondensate response to that of the condensate. 

Here, instead, we evaluate the total current response as the sum of the response of the 
noncondensate at fixed condensate and of the response of the condensate to the U field. Both 
contributions are obtained by a functional-derivative technique from the equations of motion of 
the single-particle Green's functions G . In our calculations we take the single-particle self- 
energies as functionals of G and of the order parameter <y>, this approach being the most 
general and naturally amenable to perturbative treatments [33]. The results can be immediately 
compared with those of Griffin [14], allowing us to identify the two contributions to the total 
current response as coming from the irreducible and from the reducible diagrams, respectively. 

3.1 Linear response of the noncondensate. 

We define the response function for the current density of the noncondensate through 
the functional derivative of its Green's function with respect to the potential U at fixed 
condensate (i.e. with the changes in condensate density and superfluid velocity set equal to 
zero): 

X -J\,2) = /(2m)" 2 [V u ,V 2 , 2 , f^ffyW.O W*. , (3-2) 

where V { v = V 1 — V v . The equation for the general four-point response function SG a/3 1 dU is 
obtained by taking a functional derivative on the Dyson equation for G ap : 

— ^ = G aY (l,2)GJ2,V) + G a Jl,3)— ^ -G SB (3\Y) , (3.3) 

dU(2,2') 7 7/3 7 5U(2,2') sp 

where E a/} are the single-particle self-energies, 
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z *™-8i$k?° • (3 ' 4) 

with c^l) = (7(1) and c 2 (l) = (J*(l) [34]. The convention of integration over repeated barred 
variables and of summation over repeated Greek indices has been used. Equations (3.3) and 
(3.4) are used in Appendix B to derive a microscopic expression for (see eqns (B.l) - 

(B.3)). 

Of special interest here is the relationship between and the noncondensate density- 
current response From the Dyson equation for G we find the equation of motion for the 

noncondensate density n(\) = /G u (l,l + ), and hence by taking its functional derivative with 
respect to U at fixed condensate we get the Ward identity 

X .-.(1,1' ) = {icoYW, ■ * (1,1' ) + m-W lPn (l,Y )] . (3.5) 

"J JJ 

In eqn (3.5) we have defined 

V 1 p„(l,2) = V 1 [n(l)5(l,2)]-Re{V 22 ,[< ifr\l)> eq ^^ k,} . (3.6) 

With the further definition 

p s {\,2) = n{\)8{\,2)-p n {\,2) , (3.7) 
eqn (3.6) can be rewritten in terms of the vertex functions introduced in Appendix B as 

V 1 p j (l,2) = -m^ A (l,2) (3.8) 

(see eqn (B.9)). Equations (3.6) and (3.8) extend to the inhomogeneous Bose fluid an equation 
already derived in the homogeneous case by Huang and Klein [20]. It is easily shown that the 
function #,(1,2) reduces for the homogeneous fluid in the static limit to the microscopic 

definition of the superfluid density given by Griffin [14]. Evidently, the function p„(l,2) 
reduces to the normal-fluid density. We shall make use of these functions in dealing with a two- 
fluid model in § 4 below. 

As a final remark, we notice from eqn (3.5) that in the limit indicated above p n may be 

interpreted as the so-called diamagnetic contribution to the irreducible current-current response 
function (see also the discussion given by Hohenberg and Martin [12] for the homogeneous 

Bose fluid). 
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3.2 Off-diagonal response and total current-current response. 

In section 3. 1 we have treated only the irreducible part of the current-current response. 
Detailed calculations of the off-diagonal (condensate-noncondensate) response and of the 
reducible contribution to the current-current response are given in Appendix B. The final result 
for the current-current response function Xjj can be cast in a form which is equivalent to eqn 

(6.38) in the work of Hohenberg and Martin [12]: 







® 
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A comparison of eqn (3.8) with eqn (B.15), i.e. 
* (1) I^(U) 



<5v, 

\ dtl cJ 



(3.9) 



5cp(2) 
yields the result 



V lPj (l,2) = V 2 



TO 



(3.10) 



(3.11) 



Sv t (2)° ' 

The microscopic expressions (3.9) and (3. 1 1) should be contrasted with the results 

reported by Hohenberg and Martin [12] for a two-fluid model subject to slowly varying 

perturbations, upon neglect of fluctuations in the condensate density (see their eqn (5.16)). 
These authors show that in such a limit %~ is the response of the current density j(r, t) to an 

external field given by the normal-fluid velocity \ n (r, t) and ^ v v is the response of the 

superfluid velocity v s (r,t) to an external field given by the interdiffusion current j r (r,t) = 

j(r,0-n(r)v„(r,0 =p,(r)[v,(r,f)-v II (r,0], p s (r) and p„(r) = n(r)-p s (r) being the 

equilibrium densities of the superfluid and normal-fluid components. Furthermore, their 
definition of superfluid density involves, as in eqn (3.11), the functional derivative of the 
current with respect to the superfluid velocity at fixed normal-fluid velocity. 

In summary, the main results of this section are the microscopic relations for the 
generalized normal-fluid and superfluid densities given in eqns (3.6) and (3.11), and the 
expression (3.9) for the total current-current response. These are supplemented by eqn (B.l) 
for the irreducible current-current response and eqn (B.7) for the off-diagonal condensate- 
irreducible current response. We have also made full contact with the two-fluid model as 
discussed by Hohenberg and Martin [12], for what concerns in particular the superfluid density 
and the normal-fluid velocity. With regard to the latter we also recall that Kane and Kadanoff 
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[35] have interpreted v n (r,t) as the local velocity of bodily flow of the fluid within a 
microscopic Boltzmann equation approach to the two-fluid hydrodynamic behaviour. 

4. Current functional formalism for the two-fluid model of a weakly 
inhomogeneous superfluid 

In contrast to the Hohenberg-Kohn-Sham structure that we have derived in § 2 for the 
linear response of an inhomogeneous Bose condensate, the structure of the current response 
and cross response functions given in § 3 is quite complex. However, in the last resort one 
may obtain a practicable calculational scheme from the general formalism only by (i) taking 
special limiting cases and (ii) introducing approximations. An important limiting case is that 
of weak inhomogeneity, in which both the unperturbed density profile of the fluid and the 
external fields acting on it are slowly varying in space. Recent progress for the current 
response of electronic systems in the normal state [19-23] has been referred to in § 1. 

Equations (3.5) - (3.7), introducing the functions p„(l,2) and p s (l,2) which for the 
homogeneous fluid in the static limit reduce to the densities of normal fluid and of superfluid, 
can be taken as the basis for developing a two-fluid model for the generalized hydrodynamics 
of the inhomogeneous superfluid in the linear-response, weak-inhomogeneity regime. 
Specifically, we propose that eqns (3.6) and (3.7) may provide a reasonable definition of the 
equilibrium superfluid density p s (r) in the weakly inhomogeneous case, when we take its 
static limit and the k = component of its Fourier transform with respect to rj - r, . The 
functional derivative entering the RHS of eqn (3.6) is a five-point correlation function, to be 
evaluated on the equilibrium state. 

4.1 Generalized hydrodynamics of a homogeneous superfluid. 

With this definition of the superfluid density, we first proceed to extend Landau's 
hydrodynamic theory for the homogeneous superfluid [13] to finite-frequency phenomena. 
We use for this purpose the well-known memory function formalism, as described e.g. in the 
book of Forster [28]. We assume isothermal conditions, i.e. neglect the couplings between 
density and temperature fluctuations. The form of the generalized hydrodynamic equations is 
in fact dictated by some general considerations: (i) invariance under a Galileian transformation 
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and Onsager symmetry must hold; (ii) as a consequence of the zero-force and zero-torque 
theorems, the time derivative of the current density j is driven by the divergence of a 
symmetric tensor of the second rank; (iii) the time derivative of the superfluid velocity v s is 
the gradient of a scalar quantity, in view of its irrotational character below threshold for 
vortex generation; and (iv) as already recalled at the end of § 3, the internal driving forces are 
determined by the divergence of the normal-fluid velocity v n and of the interdiffusion current 
j r = j - n\ n = p s (v s - v„) . The form that we propose for the generalized hydrodynamic 
equations thus is 

-imco8}(r, co) = -V ■ [8p(r, co)l] + V ■ <r(r, co) (4. 1 ) 

and 

-imQ)8\ s (r,co) = -V[8n(r,co)] + V-o (s \r,co) , (4.2) 

where the stress tensors are given by 

dv dv 2 r 
a, = W(0)- p (n) I ico](-^ + - - 5,V • v„ ) 

dr. 3 

+<5,[C 2 («)V-v„ + C 1 (»)V-j r ] (4.3) 

and 

C7f = 5,[C 3 («)V-j r + C 4 («)V-vJ . (4.4) 
Here, p Q (n) is the ideal-gas pressure at the equilibrium density n, while 8p(r,a>) and 
8{l(r, a>) are the local pressure and chemical potential fluctuations. These are given in terms 
of the density fluctuations 8n(r, CO) and of the entropy fluctuations &(r, co) by the linearized 
expressions 

8p(r,C0) = (nK T )~ 1 8n(r,C0) (4.5) 

and 

8/n(r,co) = (n 2 K T )- 1 8n(r,co) + (Ts I c v )8s(r,co) . (4.6) 
Here, K T , s and c v are the isothermal compressibility, the entropy and the constant- volume 
specific heat (per particle). Of course, in isothermal conditions the entropy fluctuations are 
associated with fluctuations in the density of the noncondensate. 
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The functions and r\ in eqns (4.1) - (4.4) are frequency-dependent visco-elastic 
coefficients, whose real parts express attenuation spectra. Within the memory function 
formalism these spectra are related to the response functions of the system by the following 
generalized Kubo formulae: 

4 —(dm 2 
Re[£ 2 (G>) + - 77(G))] = lim t ^ — — Im£jj(*,G>) , (4.7) 

ReT](0)) = \im k ^ Q -^-lmxl(k,0)) , (4.8) 
k 

ReC(ft)) = Um^ ^Imj ¥V (k,(0) (4.9) 

and 

Re[C 1 (ft))] = Re[C 4 (ft))] = lim^ -^Im4(fc,O)) , (4.10) 

the equality £ = £ 4 being due to Onsager symmetry. Of course, the imaginary parts of these 
spectra, which have the meaning of finite-frequency elastic moduli, are related to the real parts 
by Kramers-Kronig relations. 

The next step involves relating the visco-elastic functions to ex kernels introduced from 
the long-wavelength behaviour of the j - v s response functions and of their single-particle 
equivalents. Within the two-fluid model, i.e. setting j = p s \ s + p n \ n and n= p s + p n , we use 
the results reported in eqns (4.24) - (4.26) in the work of Hohenberg and Martin [12] and those 
given in Table 2. 1 of the book of Nozieres and Pines [36] to find the following limiting 
behaviours: 

lim k ^ Xv ^(k,CO) = (p a r8 ap +A ap (co)k 2 +o(k 2 ) (4.11) 

and 

lim^ ^ l , /j (fc,«) = (pJ- 1 <5 a/s +A a/J W 2 +o(fc 2 ) . (4.12) 
Here, the Greek indices refer to the superfluid (s) and to the normal (n) component, %l aV are 

the single-particle response functions and the A 's are frequency-dependent coefficients. The ex 
kernels are then defined through the equation 

t to 2 

f a ,(ft)) = lim^ ^p a p /J [^ ai> (fc, G ))- Zi a ^(fc, ft ))] . (4.13) 
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Use of the relations (4.7) - (4.10) and of their Hilbert transforms yields the following 
expressions for the visco-elastic coefficients in terms of the ex kernels (4.13): 

£ 1 (G);n,T) = -(i(or% LV ((Q;{pJ)-dp ex (n,T)/dn\ T ] , (4.14) 

C 2 (0);n,T) = -{icoy% Lh (co;{pJ) - Af hh {co;{p a }) 1 3 

-nd Pex (n,T)/dn\ T ] , (4.15) 
C 3 (co;n,T) = -{i(Q)- l [f v v {(Q-{p a }) - d^l ex (n,T) I dn\ T 

-(Ts/c v )dnJn,T)/dT\J , (4.16) 
C 4 (co;n,T) = -(ia>r l \f vjL (©;{/>„}) - dp ex (n,T) I dn\ T ] (4.17) 

and 

r 1 (Q);n,T) = -(iQ))- 1 f jTjT (0);{p a }) . (4.18) 
Here, p ex and fi ex are the excess pressure and chemical potential. We have denoted by 
{p a } = (p s ,p n ) the densities in the equilibrium state and left implicit the temperature dependence 
of the ex kernels. 

4.2. Local-density theory of the weakly inhomogeneous superfluid. 

The above results are easily extended to a weakly inhomogeneous superfluid in 
isothermal conditions. The general properties that we have used above to derive the structure 
of the generalized hydrodynamic equations in the homogeneous case remain valid, while a 
Ward identity is essential to relate the effect of a weak inhomogeneity on the excess kernels to 
their density dependence. This identity is obtained by modulating the densities and comparing 
the long-wavelength behaviour of the inhomogeneous response functions, to first order in the 
inhomogeneity, with those of the modulated homogeneous system [19]. 

Following the steps given in the derivation of Vignale and Kohn [20] for the electron 
fluid, we create a weak modulation of the densities of the superfluid and normal-fluid 
components given by 8p a (r) = 2£ a p a cos(q- r) at long wavelengths q— >0,with £ a « 1. 
The ex kernels in such an inhomogeneous system can be written in the form 

Q\k,r,©) = ^/k,©;{p„}) + 2f^(k + q,k,©;{p a })cos(q-r) (4.19) 
at long wavelengths k — > . The Ward identity reads 
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lim^ f^(k + q, k, m{p a }) = ^yPr ( 4 - 20 ) 

and hence eqn (4.19) yields 

q ft (k,r,©) = f^(k,«;{p a (r)}) . (4.21) 

Namely, the kernels in the weakly inhomogeneous superfluid are given by those of the 
homogeneous superfluid taken at the local equilibrium densities. 

In summary, the combination of the zero-force and zero-torque theorems with the 
Ward identity (4.21) determines the generalized hydrodynamic equations of the weakly 
inhomogeneous superfluid in the form 

-im ft)(5j(r, (O) = -n(r)V(dn(r,(0) I n 2 K T ) + V ■ a(r,(0) (4.22) 

and 

-imco8\ s (r,a) = -Wn 2 K T y l 8n(r,co)-(Ts / c v )ds(r,co)] + V- a (s \r,co) . (4.23) 
The stress tensors entering eqns (4.22) and (4.23) are given by eqns (4.3) and (4.4), where the 
visco-elastic functions are those of the homogeneous superfluid at the local equilibrium 
densities, according to eqns (4.14) - (4.18), and can be calculated from the response functions 
of the homogeneous superfluid with the help of eqns (4.7) - (4.10). The thermodynamic 
coefficients entering the ALDA terms (the first term on the RHS of eqns (4.22) and (4.23)) are 
similarly given by the corresponding quantities for the homogeneous fluid at the local 
equilibrium density, as a direct consequence of the linearization effected in eqns (4.5) and (4.6) 
for the fluctuations of pressure and chemical potential. These ALDA terms are evidently 
responsible for the first and second sound modes in the superfluid. 

5. Summary and discussion 

In conclusion it will be useful to briefly summarize the main results that we have 
achieved in this work and to indicate some directions for future work. 

The recent experiments performed on trapped vapours of alkali Bose atoms, measuring 
their collective excitations in wide ranges of frequency and temperature from the hydrodynamic 
regime to the collisionless regime, motivate the search for a unified theory of the dynamics of an 
inhomogeneous Bose superfluid. 
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In developing the basic time-dependent density functional framework for such a system 
in the linear response regime, we have focused our analysis on quantities having an immediate 
meaning, i.e. the condensate density, the superfluid velocity, and the current density. We have 
thus recognized that the essential building blocks of the response matrix of the superfluid derive 
from the two-by-two problem posed by the response of the density and phase of the condensate 
and from the irreducible part of the current response of the noncondensate. The difference 
between the condensate and the superfluid is evident from eqn (3.9), showing that the 
superfluid can be viewed as a condensate coupled to the noncondensate. 

We have used the equations of motion for the density matrix and the condensate wave 
function to microscopically define the Kohn-Sham response functions of the condensate and the 
interaction kernels entering the whole linear response matrix of the superfluid. The basic 
equations that we have derived are eqn (2.30) for the condensate response and eqn (3.9) for the 
total current-current response, combined with eqn (B. 1) for the irreducible part of the current- 
current response and eqn (B.7) for the condensate-irreducible current response. Two further 
results of our analysis are worth emphasizing: (i) the study of the dynamics of an 
inhomogeneous Bose superfluid at finite frequency requires attention to the coupling between 
the amplitude and the phase of the condensate, as related to the violation of the continuity 
equation for this component; and (ii) the microscopic, space- and time-dependent extensions of 
the density of normal fluid and superfluid emerge naturally from the theory (see eqn (3.6), (3.8) 
and (3.11)). 

The above definition of the superfluid density has allowed us to introduce generalized 
hydrodynamic equations for a weakly inhomogeneous two-fluid model in isothermal conditions 
(see eqns (4.22) and (4.23)). We have shown that in this limit the visco-elastic functions are 
related through eqns (4.14) - (4.18) to the interaction kernels of the homogeneous fluid taken at 
the local equilibrium densities. We have also shown that local-density expressions hold for the 
thermodynamic coefficients responsible for first and second sound. As remarked in § 2.2, the 
fluctuations in condensate density decay on a faster time scale than the hydrodynamic 
fluctuations and accordingly we have not included this variable in the basic set for our 
generalization of hydrodynamic theory. This point may need further study. An extension of the 
theory to include couplings between density and temperature fluctuations may also be of some 
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interest. 

With regard to approximate calculations for a Bose superfluid, the visco-elastic spectra 
can be calculated from correlation functions of the homogeneous fluid at the local equilibrium 
density through eqns (4.7) - (4.10). In particular, the sound attenuation spectrum and the shear 
viscosity spectrum in eqns (4.7) and (4.8) can be calculated in a collisionless regime by means 
of the same two-pair decoupling scheme which has been used for an electron fluid in the normal 
state [22]. The results for these spectra as functions of temperature will be published elsewhere 
[37]. One recovers in this way for the sound-wave attenuation spectrum in a dilute Bose fluid at 
zero temperature the results previously obtained by Wong and Gould [25]. One can also show 
that in this approximation the value of r]((0) is simply proportional to that of C 2 (®) [37]. 

More generally, a second-order perturbation expansion of the single-particle self- 
energies, combined with mode renormalization as in the so-called one-loop approximation 
introduced by Wong and Gould [25] provides a simple scheme allowing one to evaluate all the 
dissipation spectra in eqns (4.7) - (4.10) and the related mode-frequency shifts [38]. It turns out 
that within this approximation all these spectra can be expressed in terms of four exchange- 
correlation building blocks: a condensate kernel (the proper part of 8a ex I 8n c ), two 
noncondensate functions (the irreducible proper parts of the longitudinal and transverse current- 
current response) and a cross condensate-noncondensate term (the antisymmetric combination 
of the vertex functions A a ). Actual calculation of these four functions yields that at all 
frequencies they take the same value aside from simple multiplying factors [38]. 
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Appendix A. Relations between excess kernels of the condensate and single- 
particle self-energies. 

We outline here an alternative derivation of eqn (2.36). We start from the expressions 
(2.12) - (2.15) for the response functions in terms of the weighed Green's functions introduced 
in eqn (2.10) and use the Dyson equations for the weighed Green's functions: 

G aP = G^ + YG^®l YS ®G sp , (A.l) 

7,5 

where 

r a/j (l,l') = [n c (r 1 )n e (r 1 ,)r 1 < w'M)> eq E ap (l,Y)< ¥p (r v )> eq (A.2) 

are the weighed single-particle self-energies. 

The equations of motions for the response functions L are easily constructed from 

linear combinations of the set of equations (A.l), and the structure of eqn (2.36) is recovered 

with the following identifications: 

^ e jS(p\„=2n c (r l )n c (r r )[Si; ) -S^-^] , (A3) 

S^/SnJ^ks^ + S^-v] , (A.4) 

:(-) c(-h 



and 



where 



and 



&&„ I 8n e \, = -mte\S$ - S«] (A.5) 
Sa ex /8(p\ n =in c (r v )[S^ + S£] , (A.6) 



s™(i,v ) = -[z n (i,r)±z 22 (i,r)] (A.7) 



Sg\l,V ) = -[E 12 (U')±£ 21 (l,r)] . (A.8) 



2 

Equations (A.3) - (A.8) make explicit the connection between the excess kernels and the 
combinations of symmetrized self-energies which were introduced in the homogeneous case by 
Wong and Gould [18]. Of course, for the homogeneous Bose fluid the evaluation of the ideal- 
gas response functions and the matrix inversion involved in eqn (2.37) are easily performed. 



Appendix B. Microscopic expressions for the current response 
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We start from the definition (3.2) for the irreducible current response of the 

noncondensate at fixed condensate. From eqn (3.3) we obtain 
^ Ti a2) = j Ti (l,2) + /(2m)- 1 [V u ,G 1/ (l,3)G 51 (3',l')U 



5LJ3,3') 
x ~ 7 - T , X ,,^A\2) 



where we have introduced the notations 

-2 



^(1,2) = i(2mr [V u .V 2>2 .G 10 (l,2)G al (2 ,1' )] 1=I , 2=2 . 



and 



x <f (i,r,2) = (2iw)~ 1 [v,,. SGa ^ 1,1 \ 

^ a/3j v , j v 51/(2,2') 



(B.l) 



(B.2) 



(B.3) 



We remark that the function in eqn (B.3) reduces in the case a = f3 = 1 and 1' = 1 + to the 
noncondensate density-current response function X nj (^ ) • A Ward identity relating this latter 

function to is derived in the main text (see eqn (3.5)). 

We turn next to the off-diagonal (condensate- noncondensate) response functions. These 
can be expressed in terms of X-.~. and of the condensate response matrix introduced in § 2.2. In 



view of the symmetry property 

S<W)> \ = . SG n (2,2) 
8U(2,2) " 1 8r]\\) u 



(B.4) 



we discuss in detail only the dependence of the density and phase of the condensate on the 

potential U at constant a and # . We proceed by functional differentiation of the equation of 
motion (2.14), after supplementing its LHS by the term | d2U{\,2) < y(2) >. We define a 



vertex function A a (l,2,3) through 

*<Va(l)>, 



8U(2,2') 



-L, = 0^(1,3^(3,2,2) 



This is related to the functional derivative entering eqn (3.6) by 



A„(l,2,2)= 5(1,2X^(1,,) >„+[ 



<5C/(2,2) 



I \ J <Vc> 



(B.5) 



(B.6) 



We then find for the condensate-noncondensate current response the result 
f Sn c {\) } 

■ ^ - "^(2,3)^ 



'■ 3 ' su(3,y) 

S<p(l) 



3.3' 



8U (3,3') 



Z„ C „ C (U) Xn c9 (\,2) 

y X 9 n c (l2) * w (l,2) 



^(2,3) 



(B.7) 



/3=3' 



where 
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and 



a A (l,2) = K^)]" 1 Re[A 1 (l,2) < W'\^)>J 
^(l,2) = 2Im[A 1 (l,2)< y\r l )> eq ] , 



(B.8) 



(B.9) 



with 

A 1 (l,2) = (2m/r 1 V 22 ,[A a (l,2,2)] 2=2 , . (B.10) 
These equations are the natural extension of those given by Wong and Gould [25] for a 
homogeneous Bose fluid (see e.g. their eqns (2.26) and (2.27); see also eqn (5.8) in the book 
of Griffin [14]). 

The vertex functions introduced above in the calculation of the off-diagonal response 
also enter to determine the reducible contribution to the current response. Starting from the 
definition 



^ jj (l,2) = /(2m)- 2 [V 11 ,V 2 



SG n (XV) 



, s !») ]i=r,2=2 



2 ' 2 ' 8U{2,T) ' n 

where G ap (\,Y ) = -i< T[ yr a (l)yr}(r )] >, we have 

jW , = jW | i SG n (U') , 5<y g (3)> 
(5t/(2,2) " <5C/(2,2') K ' Vs) 8<y/ a (3)> u 8U(2,2') " 



(B.ll) 



(B.12) 



the first term on the RHS being the irreducible part which is #.-(1,2) in eqn (3.2). The 

reducible part in the RHS of eqn (B.12), from the definitions of the density and phase operators 
of the condensate in eqn (2.2), can be written in the form of a vector product, yielding 



4 red) (l,2) = (2m)- 



<5j(D 



8n c (3) u 8(p(3) 



V 5n ^ | 

2,2' r» ^ 'r? 



<5C/(2,_2' ) 
S<p(3) 



2,2' 



(B.13) 



/2=2' 



The second vector has been calculated in eqn (B.7). Using eqns (B.12) and (B.13) together 



with the symmetry property in eqn (B.4), we find the result 



(A. 




® 


Xv s v s 


7 > 


® 








\%n c v s 


%n c n c j 





(B.14) 



This is eqn (3.9) in the main text. 

Finally, using eqns (B.7), (B.13) and (B.14) we obtain the results 

TO -i„=£a2) 



8(p(2) 



(B.15) 



and 
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5j(1) -l £/ =a A (l,2) . (B.16) 



Sn c (2) 

These complete the evaluation of the current-current response. 

As a final remark, we point out that in the homogeneous limit our eqns (B.10) and 
(B.14) reduce to those found by Griffin in the dielectric formalism (see ens (8.30) - (8.32) in 
ref. [14]). This immediately allows a diagrammatic interpretation of our results. 
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